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ABSTRACT 

The nonlinear dynamics of a warped accretion disc is investigated in the important case 
of a thin Keplerian disc with negligible viscosity and self-gravity. A one-dimensional 
evolutionary equation is formally derived that describes the primary nonlinear and 
dispersive effects on propagating bending waves other than parametric instabilities. It 
has the form of a derivative nonlinear Schrodinger equation with coefficients that are 
obtained explicitly for a particular model of a disc. The properties of this equation 
are analysed in some detail and illustrative numerical solutions are presented. The 
nonlinear and dispersive effects both depend on the compressibility of the gas through 
its adiabatic index F. In the physically realistic case F < 3, nonlinearity does not lead 
to the steepening of bending waves but instead enhances their linear dispersion. In 
the opposite case F > 3, nonlinearity leads to wave steepening and solitary waves are 
supported. The effects of a small effective viscosity, which may suppress parametric in- 
stabilities, are also considered. This analysis may provide a useful point of comparison 
between theory and numerical simulations of warped accretion discs. 
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1 INTRODUCTION 

The existence of warped accretion discs, in which the orbital plane of the gas varies slowly with radius and possibly with 
time, is suggested by both observational evidence and theoretical reasoning. Of the observational results, perhaps the best 
examples are certain X-ray binar y stars, including Her X-1, in which the X-ray source appears to be periodically occulted 
by a precessing warped disc (e.g. iGerend fc BovntorJ Il97d Iciarkson et al.lf200.'y) . and a few active galactic nuclei, including 
NGC 4258, in which the warped shape of the disc is revealed by maser emission (e.g. iMivoshi et aljEoflfil : ICreenhill et al.l 
I2OO3I) . Theoretical arguments indicate that a fluid disc will be warped by an external torque if it is misaligned with the 
equatorial plan e of a spinning black hole or magnetized star at its centre JSardeen fc Pett erson[jl^7^Jljj3uii£v_fc^ 
llQSOl : lLailll99fll ). or with the orbital plane of a binary companion, planet or other satellite ijPaDaloimufc^rau^nJl^^ 
Even in an initially aligned system, th e disc may develop a warp through a linear instability of the coplanar state, depending 
on tidal, radiation or magnetic forces (lLubowlll992l: IPringljll9M: lLa]ll99fl^ . Interpreting the behaviour of warped discs is an 
important and challenging problem in astrophysical fluid dynamics. 

The basic aim of theoretical approaches in this subject has been to understand how the shape of the disc evolves under the 
action of internal stresses and external torques, and to derive practical one-dimensional equations that describe this evolution. 
When t he warp is sufficiently small t he g overning equations are lin ear and can be obtained from a perturbation analysis of a 
flat disc. lPapaloizoti fc Pringl3 l)l98,^ and^gg^lgizgj^^^in ij^^H) deri ved lineariz e d equ ations for warped discs, superseding 
the pioneering but flawed analyses of iBardeen fc Petterson il975l) and IPettersonI il97d) . and demonstrated that there are 
two basic dynamical regimes for warped Keplerian discs in linear theory. Let the disc have angular velocity 0(r) a nd vertical 
scale-height H(r), and let a denote the dimensionless effective viscosity parameter of IShakura fc SunvaevI il973l) . To a flrst 
approximation, if a>H/r the warp satisfies a diffusion-type equation with an effective diffusion coefficient H^Q/{2a), larger 
by a factor of 1/(2q^) than the effective kinematic viscosity of the disc. On the other hand, if a<_ff/r the warp satisfies a 
wave-type equation and propagates with speed Hn/2. The rapid diffusive and non-dispersive wavelike behaviours are unique 
to Keplerian discs as they result from a resonance between the orbital and epicyclic frequencies that couples the vertical 
motion associated with the warp to shearing horizontal epicyclic motions. 
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A possible limitation of these approaches is that the Eulerian perturbation analysis on which the linearized equations are 
based requires that the vertical displacement be small compared to H, and is therefore formally invalid for any observable 
warp. It is appreciated that the linear theories may be valid for larger amplitudes but this can only be demonstrated using 
a Lagrangian or semi-Lagrangian method. More recentl y it has also become possible to test these the ories and to study 
aspects of warped discs using numerical simulations fe.g. iLarwood et al]ll99fil : iNelson fc PaDalo izyiJ|l^^ ^. Such s i mulat ions 
have focused mainly on the wavelike regime a^H/r and have found reasonable agreement with IPaoaloizou fc LinI ijlQQIih for 
small-amplitude warps. There are significant discrepancies, however, which suggest that nonlinear and dispersive effects can 
be important. 

A particular concern has been that the shearing epicyclic motions associated with the warp in either regime are predicted 
to be very fast, comparable to or larger than th e sound speed, for observabl e warps. Shocks may form if the warp has a short 
enough wavelength that these motions collide ([Nelson fc Papal piz^xj [l j)9flll . but the horizon t al motions may also be para- 
metrically unstable and decay into smaller-scale inertial waves iPapaloiTOiLi fc TerauemlllQQsl : ICammie. Goodman fc Ogilvij 
I2OO0I) . Numerical simu lations indicate that th e shearing motions are also damped by magnetohydrodynamic turbulence in a 
quasi- VISCOUS manner llTorkelsson et ajboool). Further studies of all these effects are needed to assess their importance. 

In earlier wo rk (logiIv i?199^r*Ogilvic"200fl) I derived a fully nonlinear theory of warped discs with effective viscosity, 
which agrees with lPapaloizou fc IMng lg, Uggjfl in the appropriate limit and also bear s a formal resemblance to the simplified 
nonlinear equations adopted bv lPringy (|l99a ). However, the theory of IOgilvi3 (|l99flll does not describe the regime of propa- 
gating bending waves in Keplerian discs with a < H/r. It was noted there that the nonlinear dynamics of the wavelike regime 
is likely to be very different, and the purpose of the present paper is to investigate this important case. 

In Section 121 of this paper I review the linear theory of long- wavelength bending waves in Keplerian discs and introduce 
the detailed derivation that follows in Section |21 The nonlinear evolutionary equation is analysed and solved numerically in 
Section |4| and conclusions are presented in Section 



2 LEADING-ORDER DYNAMICS 

The linearized equations for long-wavelength bending waves in a thin, inviscid, Keplerian disc are (e.g. iLubow fc Qgilvielbood) 

Er^n— = -^ (1) 
dt r dr^ 

where W{r,t) and g{r,t) are complex variables describing the warp and the internal torque. Specifically W — Ix +ily encodes 
the horizontal Cartesian components of the unit vector I normal to the plane containing an annulus of the disc at radius r 
and time t, while Q — Gx + iGy does the same for the horizontal internal torque in the disc (divided by 2n). Also f2(r) is the 
angular velocity, E(r) is the surface density and H{r) is an effective density scale-height, defined by 

, E= / pdz, T.H^ ^ / pz^dz, (3) 

— oc — oc 

where M is the central mass, p is the density and z is the vertical coordinate in the unperturbed disc. For a vertically 
isothermal disc H is the usual Gaussian scale-height. 
When Q is eliminated we obtain 

^ ^GAl^_d_f 2rr2^_dW\ 
Qt2 4 J]^3/2 \ Qj. ) ' 

which implies that linear bending waves propagate non-dispersively at speed HQ/2 JPapaloizou If the model of 

the disc is such that the product E_ff is independent of r, we obtain the classical wave equation 

in the conveniently rescaled spatial coordinate 
GM^,^,Y^/\,/, , „ /■ dr 



(5) 



With plausible applications and numerical simulations in mind, we consider the case H/r = e — constant; then E cx and 



4 

— ^ (X r 

Sen 



3/2 



An outwardly propagating bending wave in this theory has the simple form 

W^ef{x-t), A=~rnf{x~t), (8) 
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where A — 2Q /(YiH^rQ) is related to the radial velocity amplitude through 

Ur = Ke(^^e-'"^Y (9) 

The function / is unconstrained except by the initial conditions, and the wave propagates without change of form. 

This long-waveleng th theory of bend i ng wa ves is subject to small corrections depending on the quantity kH, where k is 
the radial wavenumber. lLubow fc Pringlel il993l) analysed the linear adiabatic wave modes in a vertically isothermal accretion 
disc in the case \kr \ <^ 1. The bending wave of interest corresponds to m = 1 and n = in their notation, and its dispersion 
relation in a Keplerian disc can be expanded in the form 

S = - (^) + 

where F is the adiabatic index. The first term on the right-hand side corresponds to the approximation of non-dispersive 
propagation as obtained in the long-wavelength bending-wave theory, while the second term indicates the primary dispersive 
correction. The dispersion depends on the compressibility of the fluid and is a weak effect if the wavelength is long compared 
to H. 

Nonlinearity will also modify the above theory. If the nonlinearity is weak and of the same order of magnitude as the 
linear dispersion, we may expect a similar outwardly propagating wave solution to equation JHJ to exist, except that / 
will no longer be an arbitrary function determined only by the initial conditions, but will evolve according to a nonlinear 
evolutionary equation. The situation is analogous to the classical p roblem of long w ater waves, where weak nonlinearity and 
dispersion combine in the Korteweg-de Vries (KdV) equation fe.g. IWhithamlE974l) . In that case nonlinearity and dispersion 
have competing tendencies and solitary waves (specifically, solitons) are found. Another well known example is the nonlinear 
Schrodinger (NLS) equation, which arises in diverse fields including nonlinear optics. If the sign of the nonlinear term is such 
as to steepen waves, the NLS equation also admits solitons. The KdV and NLS equations are recognised as generic nonlinear 
equations that arise in many applications. 

Since the bending wave propagates through an inhomogeneous medium, it is not expected in general that its evolutionary 
equation will have constant coefficients. However, it is found below that for the self-similar disc model in which H (x r and 
E cx r~^, which has reasonable scalings and is well suited for numerical simulations, the coefficients of the equation can indeed 
be made constant by an appropriate choice of coordinates. 



3 DERIVATION OF THE EVOLUTIONARY EQUATION 
3.1 Basic equations 

The equations governing the dynamics of an ideal, compressible fluid are the equation of mass conservation. 

Dp = — pV ■ u, 

the adiabatic condition. 

Dp = -FpV • u, 

and the equation of motion. 



Du = -V<E> 



I 



Vp, 



(11) 
(12) 

(13) 



where u is the velocity, T) — dt+u-V is the Lagrangian time-derivative, p is the pressure and $ is the gravitational potential. 
For a polytropic gas F is a constant. 

Consider a non-self-gravitating, Keplerian disc around a point mass M. The time-dependent warping of the disc is bes t 
described in a nonlinear regime using the warped spherical polar coordinate system {r,6,(j)) introduced bv IOgilvi3 (|l99{i ). 
In this scheme the warped mid-plane of the disc is defined by the coordinate surface 9 — 7r/2, and the warp consists 
of a tilt angle f3{r,t) together with a twist angle "f{r,t). The Cartesian components of the unit tilt vector are given by 
I = (sin /3 cos 7, sin (5 sin 7, cos 0). 

When expressed in warped coordinates, the governing equations become 

\-^d^'V4\ , (14) 



Dp=-p 
Dp = -Fp 



-dg{vg sin 9) + ■ 



{vo sin 9) + 



1 



r sin 9 



(15) 
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Due + ^'^"^ — — [«<A cos + r (D/3) sin — r(Dj) sin /3 cos 0] = — — c^eP, (17) 

r r sin 9 pr 

Dm0 + ^^^^ H [ua cos 6 + r (D/3) sin - r(D'y) sin /3 cos 0] = ^-r—^dsp, (18) 

r rsint^ prsmt^ 

where (t;^,?;^,?;^) are the components of velocity relative to the moving coordinate system, 

V> = dt + v/dr + -de + -^^d^ (19) 
r r sm 6' 

is the Lagrangian time-derivative, 

T> = dr ~ [{drfi) cos (j) + (9r7) sin fi sin 9^ — [—{drP) cos S sin + (9r7)(cos /3 sin 6 + sin (5 cos 6* cos (p)] -^—^d^ (20) 
is a modified radial derivative, and {ur,ug,u^) are the absolute velocity components 

Ur = Vr, (21) 

ug — ve + r{Df3) cos + r^D'y) sin /3 sin 0, (22) 

Ucjy — Vcjy ~ r(D/3) cos S sin (p + r(D7)(cos /3 sin 6 + sin /3 cos 9 cos (/>). (23) 

These equations involve no further approximation and are valid for arbitrary warps. As in lOgilvi^ (Il99{i^ . meaningful dynamical 
equations for /3 and 7 can be obtained only if the disc is thin and is defined to lie close to the surface 9 — tt/2. This constraint 
is implied by the asymptotic analysis that follows. 



3.2 Asymptotic expansions 



We utilize the small parameter 5^1, such that the angular semi-thickness of the disc is H/r = e — 5^. The equations are 
to be expanded in a region of the (r, t) plane that corresponds to a neighbourhood of the line x — t followed by the nominal 
centre of an outwardly propagating bending wave at leading order. (The problem of an inwardly propagating wave can be 
considered using the time-reversal symmetry of the problem.) We therefore introduce the scaled variables ^ and r to resolve 
the vertical structure of the disc and the evolution of the wave. These are defined by 



t = 5- 



(24) 



and the equations are to be valid where ^ and r are 0(1). Thus ^ is a scaled vertical coordinate relative to the mid-plane of 
the thin disc, ^ is a scaled horizontal coordinate relative to the nominal centre x = t oi the travelling wave, and r is a scaled 
time coordinate that follows the solution over the time-scale on which the radial location of the centre of the wave changes 
by a factor of order unity. The scaling of ^ means that at any one time the wave is described within a region of radial extent 
comparable to the geometric mean of r and H (Figure 1). 
Partial derivatives transform according to 



dt = —Sde + 6 Ot, dr = 5 —de, de = S dr. 

ril 

The radius r can be expressed in terms of ^ and r and expanded in the form 



r = 2 



-4/3q2/3 



(GM) 



1/3 2/3 



2C 



3r 



1+5^ + 



9r2 



= r-o(r) + (5ri(^,r) + 5 r2(C,r) + ■ 



and this allows any function of r to be expanded similarly. In particular, the Keplerian angular velocity is 



( 



4 
Sex 



4 

37 



(f) 



+ ■ 



f7o (r) + <5fii (C, r ) + r (e, r) + 



We then propose the asymptotic expansions 
P = S^ [/3o(C,r) + 5/3i(C,r) + ■■■], 
7 = 7o(C,^) + S'yi{^,T) H , 

p = 5^^ [pe(r,C) + 5pi(0,C,C,T) + 5'p2(0,C,C,r) + ---] , 
P = 5'"+" [pc(r, C) + 5pi(0, C, C, r) + S''p2i<P, C, e, r) + ■ ■ ■] , 

Vr=S [SVrl((l)X,^,T) +S'^Vr2{cl>,(,^,T) +S^Vr3{(j>,C,^,T) + ■ 
V0 = [Sv0-i_{(j>,(,^,T) + S^Vg2{4>X,^,T) H ] , 



(25) 

(26) 

(27) 

(28) 
(29) 
(30) 
(31) 
(32) 
(33) 
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Figure 1. Typical spatiotemporal domain of an outwardly propagating bending wave. The solid line corresponds to the curve x = t for 
a disc with H/r = = 0.01. The dashed lines bound the region |§| < 10 and indicate the region that might be occupied by a wave of 
modest radial extent. Here r is expressed in arbitrary units and t in units of the orbital period at r = 1. 



v^=r{n~ D7) sin 61 + 5 [Sv^.! (0, C, ^, r) + S\^2 {4>, C, C, r) + (5'«03(0, C, C, ^) + ' ' ■] ■ (34) 
Here pc(r,() and Pg(t', C) s-re the density and pressure of the unperturbed disc in hydrostatic equihbrium, which satisfy 
9cPc = -~pcr^O.^(. (35) 
When expressed in terms of ^ and r, these quantities have expansions 

pc^ Po{C,r)+5pci{C,^,T)A , (36) 

Pc =Po(C,C,^) + '5Pci(C,C,^) + ---, (37) 
and equation H35|l is satisfied at every order in S, in particular 

dcPo = -poro^oC- (38) 



The scaling of f3 implies that the tilt angle of the disc is 0{5^) and therefore comparable to H/r, while the scaling of Vr 
implies that the radial velocity is comparable to the sound speed. It is found below that, with this choice of scalings, the effect 
of nonlinearity in the solution is of comparable magnitude to the effect of linear dispersion. The indexing of the variables in 
expansions 1281 - II34II is designed so as to give a certain order to the equations that are deduced below. The parameter a is 
unspecified and drops out of the analysis. 



3.3 Basic structure of the disc 

When these expansions are substituted into equations 114t - lll8ll . we obtain a number of equations to be considered in turn. 
Equation Ijl^jl at leading order [0(1)] is satisfied, as the Keplerian rotation balances the gravitational force. Equation at 
leading order [0(5'^)] is satisfied by virtue of the vertical hydrostatic equilibrium, equation J^HJ- 

Equation l|35|l can be solved if the relation between pressure and density is specified. In order to make analytical progress 
we consider an isothermal model for the vertical structure, such that 

Pc = Pme"'' (39) 
Po=Pme-^'''^ (40) 

where Pm{r) and Pin{r) are the (scaled) density and pressure on the mid-plane. In order to satisfy equation I35II . these 
quantities are related by 

Pm = pmr^D.^. (41) 



The surface density of the unperturbed disc is 
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E = 52^+2[Eo(e,r) + ■■■], (42) 
with 

/oo 
por-odC = (27r)'/Vopo. (43) 
- oc 

In order that S oc r~^, as adopted in Section |21 we require that pm oc r~^. 
3.4 Horizontal and vertical problems 

The equations derived at higher orders have a consistently repeating formal structure that merits analysis. The horizontal 
components of the equation of motion, equations 1161 and 11811 at 0{5"^^), n > 1, involve the unknown quantities Vm and 
v,f,n, and have the form 

Qod^Vrn — 2QoV^n = Frn, (44) 
^lod^V^n + ^QoVrn = , (45) 

where the right-hand sides involve quantities known or partially known from the lower orders. These equations may be 
combined to give 

-2noidl + l)v4,„ = Fi,„, (46) 
where Fhn is the horizontal forcing combination 

Fbn^ Frn-2d4,F4,n. (47) 

The linear operator 9| + 1, subject to periodic boundary conditions in 4>, has eigenvalues 1 — and eigenfunctions e"'""^, 
where m is any integer. The special case m = 1 is a null eigenfunction, or complementary function, corresponding to an 
epicyclic motion or eccentric distortion of the disc with an arbitrary dependence on (. The condition for equation (|46t to have 
a solution is 

e'-^Ftn dcl> = 0, (48) 



i.e. that the forcing combination should contain no m = 1 component. (Since Fhn is real there is no difference in considering 
m = 1 or m = —1.) 

The other three governing equations give rise to a vertical problem. Equation 11411 at 0(5^"^"), 11511 at 0(5^''^'*'''") and 
1171 1 at 0{5^^"), n > 1, involve the unknown quantities p„, p„ and von, and have the form 
Po 
ro 

rpo 

ro 

+ -T^^cPo = Fen. (51) 
The quantities p„ and p„ may be eliminated, using the hydrostatic condition 1381 . to obtain 

CvBn = -Fvn, (52) 

where £ is a linear operator defined by 

Cve„ = -d(:{rpod(Vg„) + porlill{dl + l)ven, (53) 
and Fvn is the vertical forcing combination 

Fv„ = rlD.l(;Fpn + rodcFpn + pQrlQ.od^F0n- (54) 

We define the eigenvalues A and eigenfunctions w of £ as the solutions of the equation 

Cw = Xpor^Q^^w, (55) 

subject to the conditions that w have periodicity 2n in (j) and be regular at the surfaces of the disc where the density and 
pressure vanish. With these boundary conditions £ is self-adjoint with weight function po- A vertically isothermal disc has 
no definite surface and the vertical boundary condition is instead that the wave energy fiux tend to zero as |("| — * oo. The 
eigenfunctions and eigenvalues are then 



O.od4,pn - 


yen „ 
ro 




Ven a 

d(Po - 

ro 




—dcPn 

poro 



= Fpn, (49) 

-dcven = Fp„, (50) 
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wi,m = e-"""^ Hef (C), \i,m = ^ r - + 1, (56) 

where I > Q and m are integers, and He denotes an Hermite polynomial.^ We will require only the first four, Heo(C) = 1, 
Hei(C) = C, He2(C) = C' " 1 and He3(C) = C' " 3C. 

The special case ioo,i = e~"^, Ao,i = is a null eigenfunction corresponding to a rigid tilt of the annulus at radius r. For 
general (irrational) F this is the only null eigenfunction. The corresponding solvability condition on equation 1521 is 

e''*Fvnd<^dC = 0. (57) 



The forcing terms required for the development of the solution to the desired order are listed in Appendix A. 



3.5 Development of the solution 

Our objective is to determine the spatiotemporal devlopment of the warp, and we must therefore obtain equations for drPo 
and 9t7o in terms of /3o, d^Po, d^^o, etc. This is done by solving the horizontal and vertical problems in turn up to the required 
order and extracting the relevant solvability conditions. Although straightforward in principle, this procedure involves arduous 
algebraic manipulations and was carried out with the aid of Mathematica. 

The first horizontal problem is unforced: Fri = -F^i = 0. The general solution is an epicyclic motion (or eccentric 
distortion) with an unknown dependence on position and time, 

«ri = (7i (C, C, ^) cos <^ + Vi (C, C, ^) sin (58) 

Based on our experience of the leading-order dynamics (Section |5J we anticipate that in the absence of viscosity the solution 
for an outwardly travelling bending wave will be Ui — — rof2o/3oC and Vi = (cf. equations|H|and|^ . This assumption simplifies 
the calculation, and is verified subsequently. Then we have 

i^ri = -ror^o/SoCcos^, w^i = irofJo,3oCsin0. (59) 

The first vertical problem involves non-trivial forcing, and we find 

^f ., = 2(F - l)(acA) cos 4> + Podao sin 0)(1 - C^) - 6/3o(95 A) cos 20 + f5odao sin2,^)C, (60) 

involving two diff'erent eigenfunctions, W2,i and ^1,2, of the operator jC. The solvability condition is satisfied, and equation 
l|52|l has the solution 



{^~T^) (^f3ocosfli-h/?oaaosin0)(l-C') + (^-^)/3o(a«/9ocos2fli-h/3o9aosin2<7i)C. (61) 



vei 
ro 

(We do not add any multiple of the complementary function Wo,i, as this would amount simply to redefining the tilt /3i.) The 
quantities pi and pi may then be obtained from equations (1491 and (1501 after an integration with respect to (p, and we obtain 

pi = pr(C,?,T) + g{i(9e/3osin,^-/3oa57ocos0) [3 - F + (F - 1)C'] C 

+ (s-Tr) ^o^^f^" '''''^^ - cos 20) (1 - C')} , (62) 

Pi = pTiC, ^) + ^ { f (9^00 sin - f3odao cos 0) [F + 1 + (F - 1)C'] C 

+ (^4r) '^''(^f^o ^^^^20 - Podao COS 20) (F - C')} , (63) 
where the axisymmetric parts are required to satisfy 

^cPi" = -pTrl^llC, - porlfloPoidao)C- (64) 

We may assume that pf^ vanishes, since it would otherwise correspond to an arbitrary redefinition of the unperturbed density 
of the disc. Then 

pT = ^^Po. (65) 
The second horizontal problem is also forced, and we find 

Fh2 



i^/3o[ae/3o(l + 3cos20)+3/3o957osin20] [F - 1 + (F -f- 1)C'] - (^i^) /3o'(95/3o cos30 + /3o957o sin30)C. (66) 



r-o^o 2F' 

^ There are two definitions of Hermite polynomials. Those used here are orthogonal with respect to the weight function exp(— C^^/2). 
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Again, the solvability condition is satisfied. (This verifies that our choice for Ui and Vi was correct; if we had not assumed 
their form in advance, the solvability condition at this order would have implied equation 1591 1 Equation 14611 then has the 
solution 

«02 = -^roPo [9«/3o (1 - cos 20) - podao sin 2<j,] [r - 1 + (F + 1)C'] ^ 



4(3 - r) 



roPo (d^Po cos 30 + Pod^jo sin 30)C 



-^(72(C,C,^)sin0+ i\/2(C,e,r)cos0. 



(67) 



The last two terms are the complementary functions, presently of unknown amplitude. The quantity Vr2 may then be obtained 
from equation 1451 . and we find 



Vr2 = 



-^roPo [d^Po sin 20 - /3o9e7o(l + cos 20)] [T - 1 + (T + 1)C'] - 2(3 1 r) ''''' ^'^ ~ ™' ^'^^^ 



1 - 



Pi 



ro{d^Po sin — Pod^jo cos (f))(^ + U2 cos + V2 sin 0. 



(68) 



.3-r, 

The forcing for the second vertical problem is very complicated, but fortunately it is not necessary to solve it in detail. 
The solvability condition alone provides an equation for the desired quantities drPo and 9t7o. This involves much algebra, 
not written out here, and simplifies to 



d^Po + i/3o9r7o - d^Pi - iPidt^jo 



— 9e / Po((72+iV2)e'^«roCdC = --J^o/3o 



2ro 



{di^Po + i/?oi957o) 



— [d^^Po - PoOaof + iPod^ao + 2i{d^Po)dao] + ^^fr^ [Pod^iPo + 2A.(9«/3o) 

1 3(r + i) 



^0 3-r 



Po{diPo)dao - P'o [%7o + i(ae7o)'] 



(69) 



Although this equation contains information about the evolution of Po and 70, it also involves the unknown quantities Pi, 71, 
U2 and V2- Further information must therefore be sought. 

Fortunately it is not necessary to solve the second vertical problem to obtain ve2, p2 and p2 in detail. The contribution 
of these quantities to the third horizontal problem, which is the last problem we consider, is of the form 



Fh3 = • ■ ■ + —dcVri - 29^ ( —dcvj,i 



(70) 

ro \ ro / 

For the solvability condition we require only the m = 1 component of this quantity. In principle both the m = and m — 2 
components of ve2 could contribute, but the m = 2 component gives exactly zero contribution. Now the axisymmetric part of 
-Fv2 is 



f:2 = -~po4^oPo [2{diPo)dao + Pod^^ao] 



2r^ + r - 2 



involving the eigenfunctions wifi and w^fi of £, and so the axisymmetric part of ve2 is 



V92 = -ttI^o [2(95/30)9^70 + Podi^ao] 

"0 



(2r2 + r-2) /r-i 



r(r + i) 



C + 



V r 



(71) 



(72) 



This information is sufficient to deduce the final solvability condition in the form 

{drPo + iPodr-10 + d^Pi + i/3i9e7o)C + —^^-—di \{U2 + iV^2) e'™l = ■ ■ ■ , 

roiZoe'T^o 

where we do not write the right-hand side in full. By multiplying this equation by poroC/So and integrating with respect to 
C, we obtain 



(73) 



drPo + iPadr'jo + de^Pi + iPid^'jo + 



Eoro^o e'^o 



1. 



i-9« / Po(f/2 + iVa) e'^" roCdC = j^o/Jo + ^{d^Po + iPodao) 



2ro 



'7r -4 

2^1 ' 

1 (i5 + 6r-r2 



{—^) [5«e/9o - /3o(9ao)' + i/3o9«7o + 2\{d^Po)dao\ + ^ (^) [Pld^iPo + (2 - r)Po{d^Pof 

P^dMdao - 7^7^^^T7^^/3o [%7o + (r + l)i(9c7o)'] , (74) 



^0 (r + i)(3-r)""^ 2^0 (r + i)(3-r)' 

which may be compared with equation 1691 . The unknown quantities Pi, 71, U2 and V2 may all be eliminated by taking the 
average of the two equations, yielding the desired evolutionary equation for Po and 70. This is expressed most compactly in 
terms of the complex tilt variable Wq = /3oe'™: 

9r 



„ 3(3r - 2) . 



ir alWor%Wo + (1 - a)l¥o%l^o* + bWaid^WoY + (1 - b)Wo\di,Wor 



(75) 
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with dimensionless parameters 

_ 2(r^ + 2r + 3) ._r + 6 
3r(r + i) ' ~3r^' 

(r-3)(r + 2) , _ 2(r - 3) 

3r(r + i) ' ^ " 3r • 

We recall that the true complex tilt variable is W{r,t) = 5^Wo{^,t) + 0{S^) 



(76) 
(77) 



4 ANALYSIS OF THE EVOLUTIONARY EQUATION 
4.1 Rescaling 

Equation 17511 is a complex nonlinear equation describing the spatiotemporal development of the warp. The first, lin ear term 
on the right-hand side corresponds to the linear dispersion evident in equation If OH derived from lLubow fc Pringlel (if 993,1 .^ 
The nonlinear terms are cubic and also depend on the compressibility of the gas. The appearance of the factor (3 — F) in the 
denominator is significant and can be t raced to the i nvolvement of the mode 101,2, which is coupled resonantly if F = 3, even 
if the disc is not vertically isothermal JOgilvieUfQQgl) . In reality F = 5/3 is the case of greatest relevance. In discussing the 
properties of equation 1751 . however, it is instructive to contrast the cases of F < 3 and F > 3. We assume throughout that 
F > f. 

There is an explicit dependence on r in equation 1751 because the bending wave experiences a varying background as it 
propagates outwards through the disc. However, this dependence can be conveniently eliminated by working with the rescaled 
time variable 

The reason that an equation with constant coefficients can be obtained is related to the self-similarity of the disc model. The 
nonlinear terms can also be rescaled by defining 

v1 1/2 

(79) 



Wo 



2(3F-2)(3-r) 



3F2 

We then have 

-i^T* = %*4-s [01*1'%*-!- (f -0)*^%** -^b**(ae4')' + {1 - b)^\d^<i/f] , (80) 
where 

s = sgn(3 - F) = ±f . (8f ) 

Equation 1801 is a derivative nonlinear Schrodinger equation (DNLS).'^ The irreducible dimensionless parameters of the 
equation, a and 6, satisfy a > 1 and 6 > f for f < F < 3, while 2/3 < a < f and f/3 < fe < f for F > 3. 



4.2 Elementary properties 

Equation 1801 is invariant under translations of T and ^ and also under the 'gauge transformation' 'i' >—t 'i'e^^, which 
corresponds to a trivial rotation of the coordinate system through an angle x about the z-axis. A further symmetry is 
^ — ^, which corresponds to a 'reflection' about the centre of the wave. The coordinate T, and the equation itself, are 
invariant under time reversal (r —t)- The trivial solution 'I' = corresponds to a flat disc. However, the equation is not 
invariant under ^' 1— > constant, which might appear to correspond to an additional trivial rigid tilt of the disc. The reason 
can be traced to equation ^ where it is seen that the amplitude Wo (or ^l/) determines the non-trivial horizontal velocities 
as well as the tilt, because the solution under consideration is a wave travelling in one direction. 



4.3 Travelling-wave solutions 

It is readily verified that solutions exist in the form of uniform travelling waves 'I' oc exp{iLoT — ifc^), where uj and k satisfy 
the nonlinear dispersion relation 

2 An additional factor of 3 appears because of the relations between t and CI and between x and r. 

^ The term DNLS refers to a wide class of equations. The present equation, in which the nonlinear terms contain second derivatives, is 
not the form of DNLS most commonly studied. 



10 G. I. Ogilvie 



tj = - (l + 2s6|*j^) fc^ (82) 
Such a wave is linearly stable to long-wavelength disturbances if 

+ [l + s(2a- 1)1*1"] > 0, (83) 

and to short-wavelength disturbances if 

4s6-f (2a-f-6- l)"i*|" > 0. (84) 

In the rescaled variables ^ and T the (negative) linear dispersion coefficient is unity. Since 2a — 1 > and 6 > in practice, 
the waves are stable if F < 3 (i.e. s = -1-1) and the nonlinear terms increase their dispersion. Therefore it may be expected 
that solitary waves are not supported, and this appears to be confirmed by the following analysis. Ifowever, if F > 3 (i.e. 
s — —1), the waves are unstable except when |*|" > l /(2a — 1) > 1, a nd the nonlinearity counteracts the linear dispersion. 

Following the usual analysis of NLS solitons fe.g. IWhithamlHoT^) . we consider more general travelling- wave solutions of 
the form 



* = F{X) exp(ia;r), X = ^ - cT, (85) 

where is a frequency and c is a wave speed to be determined. Note that c is a slow speed in rescaled coordinates and is 
relative to the basic wave speed of H^l/2. Writing F — _Rexp(i$) with R and $ real, we then find 

(1 + sR'^)R" + sRR!^ = (1 + 2shR^)R^'^ - R{c$' - lj), (86) 
R[1 + s{2a - l)R'^] +2[l + s{2a + b- l)R^] i?'<E>' = cR' . (87) 
These equations are completely integrable. Equation 18711 is linear in <!>' and has the solution 

R^<i>' = ^ + h[l + s{2a-l)R^]"' , (88) 

where h is an arbitrary constant and q — b/{2a — 1). We have 7/9 < q < 1 for 1 < F < 3, while 1 < g < 25/23 for F > 3. 
Equation 1861 may then be written in the form 

{l + sR^^JR" + sRR'^ ^ (89) 
with 

nR) = ^ + ^ [1 + - 1)R'] + ^ [1 + .(2a - 1)R^] ^"^^ - (90) 
and has the first integral 

i(l + sR^)R''^ + V{R) = E = constant. (91) 

There is an obvious mechanical analogy with a particle (albeit of variable inertia) moving in an effective potential. 

A solitary wave would have ii — > and i?' — > as X ^ ±oo. The analogous particle would roll in a potential well between 
R = and R = R+>0 such that V{R+) = V{0) and V{R) < V{0) for < R < R+. In this problem, as 7? ^ 0, the effective 
potential diverges as V ~ (c-f 2shhY j&P' R? , so this behaviour is possible only if /i = —sc/2h. In this case 

V 



mR 



^l~2[l + s{2a^l)R^Y "+ [l + s(2a- l)i?']' ^"j - (92) 

and V tends to a constant as i? ^ 0. We exclude the trivial case c = 0. It can then be shown that V is a. strictly concave 
function of R? (i.e. d?V / dlyR?)^ < 0) when F < 3 (i.e. s = -1-1 and q < 1). Therefore it is impossible for a potential well to be 
formed, and we conclude that solitary waves are not supported when F < 3. 

In contrast, when F > 3 (i.e. s — —1 and q > 1), V ha a, strictly convex function of R^ that diverges to +co as R? 
approaches l/(2a — 1), and a potential well is always formed. The squared amplitude R\ can never exceed l/(2a — 1), so it 
appears that solitary waves are supported for small amplitudes R^ < l/(2a — 1) and stable extended wavetrains for larger 
amplitudes R'^ > l/(2a — 1). An analytical approximation for the solitary waves in the case F > 3 is possible when they are 
of small amplitude 7?+ ^ 1. In this limit 

V{R)-V{0)^^R\R^ ~Rl), (93) 



with 

The solution of equation 1911 for i?" ^ 1 is then 
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R « R+sech{kX), 



k = - 




(95) 



and furthermore $' ~ c/2. For a specified height (-R+) and width (1/fe), the frequency lj and wave speed c are determined by 
these relations. 

4.4 Numerical solutions 

Equation 1801 was solved numerically by applying a spatial discretization on a uniform grid, based on second-order centred 
differences, and solving the resulting coupled ordinary differential equations in time using a fifth-order Runge-Kutta method 
with adaptive stepsize. Illustrative solutions for the case F — 5/3 are shown in Fig. 2. The initial profile is a real Gaussian of 
relatively large amplitude. The solution of the linear Schrodinger equation, obtained when the nonlinear terms are neglected, 
can be written in closed form and is shown in the right-hand panels; it exhibits some linear dispersion. When the nonlinear 
terms are included, as shown in the left-hand panels, the wave broadens much more rapidly and there is a greater emission 
of short waves. (If their wavelength is comparable to or shorter than H it is likely that they will behave differently from 
the predictions of equation 1801 . The computational domain is large enough that the boundary conditions do not affect the 
solution significantly for the ranges of ^ and T shown. It was also confirmed numerically that small-amplitude solitary waves 
in the case F > 3 propagate according to equation 19511 . 

4.5 Inclusion of a small effective viscosity 

It is straightforward to generalize the analysis to include a small effective viscosity. Suppose that the effective shear viscosity 
is given by = ap/Sl, where a is independent of z and comparable to H/r. Let a = S^ag with ao = 0{1). Then the viscosity 
appears in the third-order horizontal forcing terms Frs and F^s, where its effect is to replace dr by {dr + ao^lo). Indeed, any 
process that damps the shearing epicyclic motions at a rate aSl would have an equivalent effect. An effective bulk viscosity of 
comparable magnitude would not affect the analysis. 
The final equation for Wo is modified to 



When <& is small the effect of viscosity is to cause the solution to decay by a factor T relative to the inviscid solution. 

In general, however, the nonlinearity of the equation means that it is not possible to find such an integrating factor. 



5 CONCLUSIONS 

We have derived a one-dimensional equation describing the evolution of weakly nonlinear and dispersive bending waves in a 
thin Keplerian accretion disc with negligible viscosity and self-gravity. The nonlinear and dispersive effects both depend on 
the compressibility of the gas through its adiabatic exponent F. Perhaps surprisingly, in the physically realistic case F < 3, 
nonlinearity tends to broaden the waves and enhance their linear dispersion. This is in contrast with the (hypothetical) case 
F > 3 in which nonlinearity counteracts linear dispersion and solitary waves are supported. 

Linear theories of warped accretion discs iPaoaloizou fc Pringyfl98i:IPapal OIZOU fc LiJl995ll previously made clear the 
distinction between the diffusive (a>///r) and wavelike («< Jf/r) regimes in Keplerian discs. Since they were based on an 
Eulerian perturbation analysis they could not predict the amplitude at which they would break down. In this paper it is found 
that nonlinear effects cause a significant modification of the bending wave (as the radial location of the wave changes by a 
factor of order unity) when the dimensionless amplitude of the warp |dVK/dlnr| is O^H/r)^'"^ and the horizontal velocities 
are comparable to the sound speed. 

Unfortunately it is difficult to explain in simple terms how these critical scalings arise or why the nonlinearity is steepening 
only when F > 3. The nonlinearity arises through a mode-coupling process and is complicated further by the resonant nature 
of the epicyclic oscillations in a Keplerian disc, arising from the coincidence of the orbital and epicyclic frequencies. 

The analysis in this paper does not by any means present a complete description of the nonlinear dynamics of the 
wavelike regime. We have focused on the evolution of a travelling bending wave and have not considered the interaction 
betwe en inward and outward wav es, nor the important case of steady warps maintained by an external torque. Further- 
more, iNelson fc Papaloizoul lll999l) noted that shocks can form if the wavelength of the warp is short enough that the 
shearing epicyclic motions from different parts of the wave collide. This effect typically requires |dVF/dlnr| to be 0(1). 



drWo ^ ■ ■ ■ - -^Wo, 
3r 

where the dots indicate the inviscid dynamical terms found previously. The DNLS equation 18011 is modified to 



(96) 




(97) 
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Figure 2. Numerical solution of equation 1801 in the case V = 5/3 starting from a real Gaussian initial condition. Left: solution including 
the nonlinear terms. Right: solution excluding the nonlinear terms. The real and imaginary parts of ^ are shown as solid and dotted 
lines, respectively, at times (from top to bottom) T = 0, 1, 2 and 4. The computational domain is much larger than shown here. 
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iGammie. Goodman fc Qgilvia (1200(3) also noted that the warp could decay through a parametric instability. The critical 
amplitude for this effect is proportional to the effective viscosity of the disc. 

It would be interesting to compare the predictions of this analysis with the results of numerical simulations of warped 
discs. If parametric instability sets in there should be a noticeable deviation from equation (1801 . The best way to make this 
comparison would be to set up a vertically isothermal disc with H oc r and E oc , to introduce a warp and the accompanying 
horizontal velocities l|59^ , and to follow the evolution of the complex tilt W{r, t) . The comparison with solutions of equation IjHO^ 
could then be made by transforming variables from W(r-,t) to ^'(^, T), or vice versa. 
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APPENDIX A: FORCING TERMS 

For reference, we list here the forcing terms required for the development of the solution to the desired order. 
2po 



pi 



(Al) 



p2 



-drpO 



rivgi 



(l-^)9,(Pc. + M-b + (l-^)% 

'iP-d(;po + (Pcl + Pi) I T^l^et^rl - (957o)9,^,t;rl] + —d(;Vei\ 

L ''o"o ^■o J 



31 + —dc{pcl + pi) 

ro 



+P< 



roQo 
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[a,..,-(9,7o)a,...]}, 



Fp2 = -drPO + (l - ^) ae(Pel - fil + (l 



2Vrl 

rofio 



31 + — (?c(Pel 



_rr!^^ +r(pci +pi) ( T^fc's'Vi - {dt^-yo)d^Vri] + —di;Vei\ 

{2 
roS2o 



Frl = 0, 



Vrl + 2QlV^l 



{d^Po cos (p + Pod^'yo sin (f>)d,;po, 



Frz = -a 



2vr2 2vrir-i 21)^1^1 \ ?Jei 

^T; 7W "€"''1 ^ "i'"^'^ 



, , ve2 vexTi , 



02 + ^ - a,7o + f 1 



2«ri \ „ / 2vr2 2t;r.iri 211^1^1 , „ 

c'e7i - ( ZZ7^ - 177^ I ^eTo 



rofio 



yrof2o rgSlo rof^S 



+2f2iv^2 + 2f22W<#.i + — ro^oC^ 
ro 

2 

7^ {9j (pei + Pi) + {d^Po COS + /3o5£7o sin (f))dc{pei + Pi) 

poro"o 

+[5{/3icos0+ (/3o9f7i +/3i9f7o)sin0]9fpo - (5£7o)a,^pi} 

Pel + Pi 



+2 



H 1 (9e/?ocos0 + /3o9j7osin0)a^po, 



Porofio PoTq^o porofll ^ 
Fgi = 2(2t;ri — rofio)(5{/Jo sin(/> - PodiJo cos (j>) — 2{dci>Vri){d^Po cos<p + Pod^'jo sin< 

2Vr 



Fe2 = (1 ~ rT7"l ^S""! + "^C^"! 



rono 



ro 



- + ^^i^ ) dcpi + -^dcpo + 2QovM 

\por^ Poro J PqTo 

+2(2wri - roQ,o)[d(i^l3i sm(j) - {Podi^-yi + Pid^jo) cos0] 
+2(2wr2 - ro^^i - ri^o) (9^/30 sin (/i - /3o9e7() cos<t>) 

—2{d^Vri)[d^l3i COS0 + {Pod^ji + Pid^jo) sintp] + 2TxMo{dT(3o sin0 — /3o9t7OCOS0) 

-7^ |fio9rf,Wr2 - f 1 - —TT-') [dfVri - (5f7o)90Vri] - —drVri\ (9* /3o COS (/> + /3o5«7o sin 0) 
ilo I \ roilo/ ro J 



-ro 1 



(l - ^) {[d((l3o - f3o{daof] COS <^ + [/3oa«57o + 2(a«/3o)a57o] sin 4>}, 



F^i = 0, 
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2 2ro 
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